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ABSTRACT 



I 

We extend the analysis of iPenev et al. 



(120071 ) to calculate effective viscosities 



for the surface convective zones of three main sequence stars of O.775M0, O.SSM© 
and the present day Sun. In addition we also pay careful attention to all normal- 
ization factors and assumptions in order to derive actual numerical prescriptions 
for the effective viscosity as a function of the period and direction of the external 
shear. Our results are applicable for periods that are too long to correspond to 
eddies that fall within the inertial subrange of Kolmogorov scaling, but no larger 
than the convective turnover time, when the assumptions of the calculation break 
down. 

We find a significantly anisotropic viscosity, scaling linearly with the period of 
the external perturbation and magnitudes of the differ ent components between 



three and ten times smaller than the 



ZahnI (11966 



19891 ) prescription. 



Introduction 



Turbulent (eddy) viscosity is often considered to be the main mechanism responsible 



for dissipation o: 



' tides and oscillations in convection zones of cool stars and planets 



( IGoodman fc OhI (119971 ) (from now on GO), and references therein). Currently existing 



descriptions have been used, wit h varying success , to explain circularization cut-of f periods 



for main sequence binary stars (jZahn fc Bouchet 



red edge of the Cepheid in stability strip (iGoncz 



( iGoldreich fc Keeley 



1989 



Meibom fc Mathieu 



20051), the 



19821 ) and damping of solar oscillations 



19771 ). However, this hypothesis has been far more successful in 



damping oscillations than damping ti des, and dif fere nt mechanisms 

3 ( 2005J bl):[ Ogilvie fc LinI 



the latter, especially for planets (see 



lave been proposed for 



(I2OO4J ) and references 
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therein). 



The standard treatment is to assume a Kolmogorov spectrum in the convection zone 
and apply some prescription to model the effectiveness of eddies in dissipating the given 
perturbation. Two prescriptions have been proposed to describe the efficiency of eddies in 
dissipating perturbations with periods (T) smaller than the eddy turnover time (r). 



The first prescription, due to 



ZahnI (11966 



19891 ). assumes that it is always the largest 



eddies that dominate the dissipation and that they lose efficiency linearly with decreasing 
period: 



mm 



Where i^max is some constant which depends on the mixing length parameter and is 
the local convective turnover time (or the turnover time of the largest local eddies). This 
prescriptio n has been tested against tidal circularization times for binaries containing a 



giant star (jVerbunt fc Phinney 



19951 ). and is in general agreement with observations. 



he second prescription, due to lGoldreich fc NicholsonI (119771 ) and iGoldreich fc Keeley 



(Il977l ). assumes that eddies with periods longer than T/2ii do not contribute to the 



dissipation. In that case, for Kolmogorov scaling: 



V = Vrnax mm 



T 

2'KTr 



(2) 



This prescription has bee n used successfu l 



Goldreich fc Kuman (119881 ) and 



Goldreich et al. 



V bv 



Goldreich fc KeelevI (Il977h : 



(119941 ) to develop a theory for the damping 



of the solar p-modes. If the more effective dissipation was applied instead, dramatic 
changes would be required in the excitation mechanism in order to explain the observed 



p mode amplitudes. However, this i nefficient dissipation is inc onsistent with o 



3served 



tidal circularization for binary stars (IMeibom fc Mathieu 



20051). Additionally, 



Gonczil 
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( 119821 ) argues that for pulsating stars the location of the red edge of the instability strip is 
more consistent with Zahn's description of eddy viscosity than with that of Goldreich and 
collaborators. 



GO gave a consistent hydrostatic derivation of the convective viscosity, using a 
perturbation approach. For a Kolmogorov scaling they obtained a result that is closer to the 
less efficient Goldreich & Nicholson viscosity than it is to Zahn's, providing a more sound 
theoretical basis for the former scaling. Of course, the observational problem of insuffi- 
cient tidal dissipation for stellar pulsations and binaries remains unresolved, as GO point out. 



Both 2D and 3D numerical simulations of the solar convection zone have revealed 
that t 



1989 



le picture of a 



Robinson et al. 



kolmogorov spectrum of eddies is too simplified ( IStein &: Nordlund 



20031 ). The simulations showed that convection proceeds in a 



rather different, highly asymmetric fashion. This suggests that the problem of insufficient 
dissipation may be resolved by replacing the assumption of Kolmogorov turbulence with 
the velocity field produced from numerical simulations. More importantly, an asymmetric 
and non-Kolmogorov turbulence might dissipate different perturbations differently, i.e. 
depending both on the frequency and geometry of the perturbation. S uch simulations hav e 



been used to develop a better model for the excitation of solar p-modes (ISamadi et al 



20031). 



In 



Penev et al. 



(120071 ) we reconsidered the problem of tidal dissipation in stellar 
convection zones of solar-type stars by applying the approximation developed in GO to the 
turbulent velocity field from a realistic 3D solar simulation and showed that, the scaling 
predicted by this procedure is very close to linear. The shallower scaling is explained by 
the fact that on the time scales captured by the simulation the largest eddies have typical 
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sizes comparable to the local pressure scale height and, in this regime, the velocity power 
spectrum is much shallower than Kolmogorov. 

In this paper we apply a more complete version of the same scheme to three stellar 
convection models appropriate for stars with masses of O.TTSM©, O.SSM© and IMq. We 
also pay significantly more attention to the normalization and the approximations which we 
introduce in addition to GO. 

2. Method 

2.1. The Perturbative Expansion for Discretely Sampled Velocity Field 

We follow the procedure outlined in GO and assume an external perturbing velocity 
field given by GO equation (8): 

Vt = A{t) ■ X, (3) 

where A{t) is some matrix that depends only on time, and not space. This is appropriate 
when the spatial dependence of the perturbation is on scales much larger than our simulation 
box (e.g. tides). The matrix A{t) is assumed symmetric, since the antisymmetric part 
corresponds to rotation, and is not expected to contribute to the energy dissipation. 

Introducing this velocity field will modify the convective flow (vq). The time evolution 
of the change in the turbulent velocity (Sv) due to the presence of the above external field 
can be written as in GO, equation (19): 

dtSv{x, t) — —2A{t) ■ Vo(x, t) —Vo • V(5v — 6v • Vvq — (pressure term). (4) 

GO assumed incompressibility, and hence the pressure term simply maintains that 
V • 5v = 0. We use the output of fully compressible simulations, so for us the pressure 
term should be much more complex. However, the only type of compressibility that we 
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can reasonably incorporate in the analysis is that due to the stratification of the convective 
layer, and even that we approximate by assuming a constant density scale height. This is 
discussed in more detail in section 12.31 

We then follow GO in writing equation (jl]) in Fourier space. However, since we are 
dealing with discretely sampled data, we use discrete Fourier transforms: 

1 V — -\ — o„,Y >•!■ I A"n . vn . <!>p\ 

Av, = > ,p^''\N^^ Ny^ N;,^ NtJ 

V - ^ V 9 ^ + ^ + ^ + (^) 



1 



A{t) = - ^A{n)e-'''' + A{-n)e'' 

where 27t/Q is the period of the external forcing. 

For more details on how exactly the Fourier transform is applied in the radial and time 
directions see section. 12.41 

In Fourier space to first order in A (the strength of the perturbation) and Qt^ (the 
ratio of perturbation time scale to convective turnover time) equation (jlj) is written as: 



^Vx,ti,u,(f) — ^X,ii,u 

^4> 



A(^]) • v(cj^ - fi, kA,^,^) + A(-fi) • v(cj^ + fi, kA,^,^) , (6) 



where Px,^,v = I - \i! \^^^i,'k! \^^^i,/k'^^^^^, with k'A_^,^ = k\^^^u + iz/Hp, is the discrete version 
of the projection operator GO define that imposes compressibility due only to a constant 
density scale height. 

We can then express the average rate of work done (per unit volume) on the turbulent 
velocities by the tide to lowest non-zero order as: 

^ T ^ ^ * 1 o 2* 

' ~ N"^N ^ Pu~u''^x,fi,u,p^x,fiyVx^fj,y^p' 

1 (7) 

£{n = 2itR/T) = Re {Sr^^r + Sr^r} + > Im {Sr+R,r-R + Sr+R,r+R}, 

^ — ' vrr 
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where M = N^NyNzNt. In keeping with GO, we have assumed that all frequencies have an 
infinitesimal imaginary part, which gives rise to the first term above. The second term is 
entirely due to the density stratification in the box: in the case of p{z) = const it is zero. 
This term is the most important difference between this calculation and GO. 

In deriving equation ([7]) we have assumed that the density is only a function of depth. 
Keeping the radial dependence is necessary because the simulation box encompasses several 
density scale heights, while the horizontal and temporal dependence of the density is a 
much smaller effect, entirely due to the turbulent fluctuations in the box and if we could 
average over different realizations of the turbulence they would not be present. 



2.2. Anisotropic Viscosity 

In order to extract an effective viscosity, we need to express the energy dissipation rate 
that would occur in the presence of actual anisotropic viscosity. 

Most generally, viscosity is a fourth order tensor relating the strain, given by A(t) in 
this case, and the viscous stress: 

With this deflnition the time averaged dissipated power is given by: 

£visc{n) = \ r dzKijmn{z)Re {Aij{n) A*^,^{n)} , (9) 
^ Jo 

where, in order to remain consistent with the Fourier transform conventions we simply 
replace the integral with a sum. To get the different components of Kij^n we evaluate 
equation ([7]) with A having nonzero elements at different locations, and use the above 
equation to flnd the respective viscosity coefficients. 
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The viscosity tensor {Kij^n) obeys a set of symmetries that dramatically reduce the 
number of independent components. Since the strain rate is symmetric by definition, and 
the stress must be symmetric in order to keep the viscous torque on infinitesimal fluid 
elements finite, we must have: 

■^ijmn — -Kjimn — -^ijnm- (10) 

In addition, the only distinct direction in the problem is that of gravity (i), so we expect 
the viscosity tensor to be symmetric with respect to rotation around the vertical axis. 

With all these symmetries we are left with only six independent components of Kij^n' 
-f^iiii, -f^i2i2, -f^i3i3, -f^ii33 and i^^33ii. Siucc the last two of these always appear 

together in the expression for the energy dissipation we will assume them to be equal. The 
remaining non-zero components can be found from those as follows: 

-f^2222=-f^llll5 
7^1122 = -f^2211=-f^llll ~ 2/^1212, 
^1221 = -f^2121 = -f^2112=-f^l212, 
-^^3131 — — i^l331 — K2323 \ ^^^"^ 

-^^3232 = -^3223 = -^2332 J 

-f^2233=-f^ll33) 
-f^3322=-f^3311- 

A more physically meaningful set of five viscosity components can be found by noting 
that under these symmetries the strain rate has only four distinct components: 

A = ^i + ^22; A'^^3; A = ^3 + ^^3; A = Ai - ^22 + 2iAi2, (12) 

along with their complex conjugates A-m — which transform under rotation by angle 
around the z axis as Am e^^^^Am- 
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Clearly then, if Smsd^) is to be invariant under such rotations, it must be of the form: 



1 /"^ 

Svisci^) = - dz [4:Ki\Ai\'^ + K2\A2\'^ + KqA^ + i^o'^o' + SiTocAA'] , 
^ Jo 

where the five new viscosity coefficients can be expressed in terms of Kij^n as follows 



Kq = - (iTnii + K1122) 



(13) 



K3333 
K1133 



(14) 



Kl — -^1313 

K2 = - (-^1111 — -^1122) 

Since in equation iQ we allow the viscosity to depend on depth, and there is no way 
to constrain this dependence, we have to choose some radial profile a priori. Our choice is 
motivated by mixing length theory: 



(15) 



where are dimensionless constants, that depend on the frequency of the external shear 
(Q), and Hp is the local pressure scale height. This is reasonable, since the turbulent 
viscosity should scale as some length scale times some velocity scale. Clearly the relevant 
velocity scale is that of convection, and in accordance with mixing length theory, we use 
the mixing length as the length scale, which is assumed proportional to the pressure scale 
height. If the mixing length is really the relevant quantity, we expect that the value of 



will be proportional to the mixing length parameter for the particular simulation. 



same scaling has been assumec 



1989 



Goldreich fc Keeley 



1977 



br all previous effective v i scosity prescriptions 



Goldreich fc Kumar 



1988 



Goldreich et al. 



199J). 



Zahn 



his 



1966 
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2.3. The Pressure Term 

In deriving the expression for £ (Eq. [7]) we assumed that the perturbation to the 
convective velocity field due to the tide will be anelastic: V ■ 5v + Vz/Hp = 0, with 
Hp = const. In this section we define two diagnostics which measure how important the 
ignored compressibility is. 

There are two sources of compressibility in the convective fiow: 



1. The convective fiow carrying parcels of matter through layers of different hydrostatic 
pressure, or in other words due to the stratification. 

2. Localized compression due to a possibly supersonic fiow, e.g. shocks. 

Ignoring the second one is justified, as long as the fiow velocity is much less than the local 
speed of sound. In the simulations we use, that condition is met by the unperturbed fiow 
for most of the box, with the exception of the supersonic driving region near the top. If the 
unperturbed fiow is subsonic and hence incompressible, the perturbations due to a "small" 
external field can safely be assumed incompressible as well. To measure the compressibility 
in the simulation box we introduce the parameter: 



V + z- 



•vo, (16) 



dz 

Where, p is the density averaged over horizontal slices and time. 

This quantity deviates from zero due to localized, transient compressions (e.g. shocks). 
Since those are unlikely to live longer than a convective turnover time, this quantity is a 
suitable diagnostic for the importance of such effects. 

Because we are measuring the mass averaged dissipation, in order for the perturbative 
treatment discussed above to be valid, we can only have ^ ^ 1 for a negligible fraction of 
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the mass. In figure [T] we plot the time averaged fraction of mass that resides in regions 
of the convective box which have C, greater than some value. The value of the convective 
turnover time, necessary to evaluate ^, was calculated as Tc = FWHM(f^)/ m8ix{vz), where 
FWB.M{vz) is the thickness of the layer over which Vz > max(f^)/2. 

As we can see, in all cases, ^ > 1 for less than 1% of the mass. This compression is 
concentrated near the top of the box, where the density, and hence the dynamic viscosity, 
is small. So, even though compressibility and shocks are important in determining the 
fiow that develops, no appreciable dissipation occurs in strongly compressible regions. The 
situation is further improved by the fact that we apply a window in the vertical direction 
that significantly reduces the importance, and completely ignores part of the compressible 
driving region near the top of the box in determining the effective viscosity (see section [2^ . 

We partially treat the first source of compressibility discussed above by, imposing 
Hp = const in the continuity equation. However, this is not valid for most astrophysically 
interesting convective zones. In fact in all cases considered here the density scale height 
varies by a factor of a few between the top and the bottom of the convective layer. In 
some sense, assuming Hp = const is not any better than assuming Hp = oo. We argue that 
ignoring the stratification from the continuity equation is a reasonable approximation. 

The simplest way to justify this is to repeat the evaluation of the viscosity with 
different values of Hp within the range encountered in the convective layer of interest. 

We can also gauge the importance of the stratification by comparing dlnp/dz to 
d\ia6vz/dz. We evaluate dlnp/dz directly, and estimate: 



The last expression comes from equation [6l and is correct when the last row of A(t) 
contains only a single non-zero entry. Since those are the only cases we use, this expression 




(17) 
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Mass Fraction 

Fig. 1. — The compressibility of the unperturbed flow as a function of depth for the three 
simulation boxes (O.775M0 - top left; O.85M0 - top right; I.OMq - bottom). The horizontal 
axis gives the time averaged fraction of mass with compressibility greater than the vertical 
value. 



is sufficient for us. 

In figure [2] we compare dlnp/dz to d\nSvz/dz (estimated as in the above expressions). 
We see that the logarithmic gradient of the density is approximately two orders of magnitude 
smaller than the typical logarithmic velocity gradient and hence we are justified in ignoring 
it. 
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le-04 0.001 0.01 0.1 1 10 



dlnp/dz < dv^/dz > / <Vz> 

Fig. 2. — The logarithmic gradient of the density compared to the logarithmic gradient of 
Vz, estimated as explained in the text for the three simulation boxes we considered: O.775M0 
- top left; O.85M0 - top right; LOM© - bottom. 
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2.4. Radial and Time Fourier Transforms 

Using discrete Fourier transforms to represent a data set, forces tlie assumption tliat 
the data is periodic in all dimensions. While this holds for each horizontal slice, it is violated 
for the radial and time dimensions. Ignoring the problem leads to artificially introducing 
spectral power at the highest frequencies because of the jumps at the boundaries. To avoid 
this, we need a special way to deal with the non-periodic directions. 

The usual solution is to window the data so that it goes smoothly to zero at the edges 
of the domain. This has the effect that it makes the values near the center of the domain 
relatively more important than those near the boundaries. Incidentally, this is exactly what 
we would like in the radial direction, since the flow near the top and bottom is affected by 
the artificial boundary conditions and is not representative of the actual fiow that would 
occur in a star. 

Further, as discussed in section 12.31 we expect that the compressibility of the fiow that 
we neglect might be significant in the upper end of the box, where the density is small 
and the fiow is supersonic. So making this region's contribution to the overall dissipation 
small is exactly what we would like. In fact, in the radial direction we go a step further 
and limit the window to completely exclude some part of the box near the top and bottom 
boundaries (see equations (ITSl) and (fT9!) ). 

In the time direction, as long as the time interval we have simulated is "representative" 
of the actual convection that occurs in a star, weighing the center of the interval more than 
the edges should not be a problem. 

To confirm that the chosen window function is not affecting our final result, we derive 
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the effective viscosity coefficient using two common windows: 

2 



Welch : M 
Bartlett : M 



1 - 



2t-T 
T 

2t-T 



max 



0,1 



2z - L, 



T 



max 



0,1 



aL 
2z-L 



(19) 



where A/" is a normahzation factor numerically equal to the inverse of the average of the 
squares of the window function at all the grid points, and a is a parameter determining 
what fraction of the radial span of the box we include in the analysis, that is we exclude 
(1 — a) fraction of the linear size of the box, half from the top and half from the bottom. 



3. The Stellar Models 



The three models used, represent the top 7-9 pressure scale heights of the convective 
zones of the present Sun, a 0.775 Mq and a 0.85 Mq stars. Table [1] shows the position of 
each model in the log g — log Tes plane. The full details of the n u meric al scheme and the 



properties of the solar simulation are discussed in 



Robinson et al. 



(120031 ). For a comparison 



between the mode 



Kupkal fl2005h and 



s used here and the work o f other groups, as well as observations see 



Hillebrandt fc Kupkal (120091 ). Here we present very briefly only the most 



important aspects of the models. 

The simulation boxes have periodic side walls and impenetrable top and bottom 
surfaces with a constant energy flux fed into the base and a perfectly conducting top 
boundary. The imposed flux was computed from a corresponding ID stellar model with 
the chosen mass and age, thus was not arbitrary, but the correct amount of energy flux 
the computational domain should transport outward in the particular star. The initial 
conditions of the 3D simulations were also derived from the same ID stellar models used to 
calculate the required flux. 
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3.1. Starting Models and Input Physics 



The ID stellar m odels used to initializ e each run were computed with the YREC 



stellar evolution code fiGuenther et al 



19921 ). They were calibrated to the Sun and evolved 



from the ZAMS. Both t he ID and 3D codes use the same realistic physi c s as d escribed by 



Guenther &: Demarqud (Il997l ). most notably the 



Alexander fc FergusonI (119941) opacities 



at low temperatures, the OPAL opacities and equation of state ( llglesias fc Rogers 



19961 1. 



hydrogen and helium ionization and helium and heavy element diffusion. 



Some details of the three models are given in Table [H The fractional radius is given 
as R/R0, where R is radius of the stellar body and Rq is the radius of the Sun. Both are 
defined at the point where T = Tgfj. The surface gravity and effective temperature are in 
c.g.s. units. 



3.2. Box dimensions 

The horizontal dimensions of each computational box (column 5 in Table 1) were 
estimated by assuming that the granule size will scale roughly inversely with g. The final 
column gives the number of grid points in the two horizontal and vertical directions in the 
square based box. 



4. Results 

As discussed in section 12.21 the anisotropic viscosity can be parametrized by five 
independent components. Assuming equation (IT3i) . we evaluate those components using the 
two window functions of equations (ITSll and (fT9l) each with two different values of a: 0.8 
and 0.9. In addition we use two values for the density scale height in each case: Hp = 00 
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and the volume average density scale height. 

The reason for using the volume averaged value of Hp instead of the mass averaged 
is that, this way, relatively more weight is given to the less dense top regions where the 
density scale height is small, resulting in a larger range between the two cases we consider. 

The frequency dependence of the five viscosity coefficients (ii'o, Kqi, -ft'oo'; -^i; and K2) 
is presented in figure [3] the curves correspond to a Welch window (Eq. [T8l) with a = 0.9 and 
volume averaged density scale height and the error bars show the span among all the cases 
for which we evaluated the effective viscosity. 

The fact that the error bars in figure [3] are small shows that indeed the choice of the 
window function is not important and that ignoring the depth dependence of the density 
scale height, and in fact the stratification altogether in the continuity equation is a valid 
approximation. 

We see that the same qualitative characteristics hold for the estimated dissipation in 
all 3 of our simulation boxes: the Kq> component is always approximately four times larger 
than the Ko, Ki and K2 components, which are in turn roughly four to five times larger 
than the f^oo' component and the scaling is approximately the same for all components, 
close to the linear scaling proposed by Zahn. 

Quantitatively, the effective viscosity we calculate can be written as: 




(20) 
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1 2 5 10 20 1 2 5 10 20 




1 2 5 10 20 

Fig. 3. — The five viscosity components, scaled by {v^)^^'^ Hp for tlie tliree stellar models 
(O.775M0 - top left; O.85M0 - top right; IMq - bottom). The lines correspond to the average 
of all curves representing the same component calculated using a = 0.9, Welch window (Eq. 
[T8l) and the volume averaged density scale height. The error bars correspond to the spread 
found among all the curves corresponding to the viscosity component with different windows, 
values of alpha and density scale heights. 
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where the parameters /Cm and A take the values: 



A: 


=1.2 ±0.3, 




=0.022 ±0.003, 




=0.080 ±0.01, 


/Coo'- 


=0.0046 ±0.000. 


/Cr 


=0.024 ±0.003, 


IC2 


=0.018 ±0.003, 



with the errors corresponding to the range of values encountered for different windows, 
values of a, density scale heights and stellar models. The above values were derived by 
performing a least squares fit of equation [20] to the calculated curves for T < O.Stc. The 
reason for restricting the fit to short periods is that at long periods we do not expect the 
perturbative calculation used in this work to be applicable and hence the derived slope is 
an artifact of the model rather than having any physical significance. 

Also there is no appreciable difference between the models of the different stars. The 
spread in the dimensionless effective viscosities for the three models is not much larger than 
the error bars at all frequencies, except the high end tails, where the effects of the finite 
resolution and time sampling become important. This suggests that at least for the range 
of conditions encountered in the convective zones of low mass stars the dissipation efficiency 
is not strongly dependent on the details of the convective flow. 



4.1. Anisotropy 

The above splitting of the viscosity in five components was done in order to allow for 
anisotropic dissipation. It is interesting to see how anisotropic the derived effective viscosity 
really is. The general isotropic case has only two viscosity components: a bulk viscosity {Q 
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and a shear viscosity (c.f. ?). In terms of those the isotropic Kijmn tensor is: 



From this it can be seen that the five components of the viscosity we calculate must obey 
the relations: 

Ki = K2 = v 

Ko = Koo' + (23) 

We see that the first two of these are clearly satisfied by the viscosity coefficients of equation 
[2T] to within the quoted uncertainties. The degree to which the last equation is not satisfied 
is: 

- 1 = 0.74 ±0.25. (24) 



Considering the fact that the fiows in our simulation boxes are not exactly like 
those inside stars, and the loosely estimated errors in equation [21] we can conclude that 
the effective viscosity we find is only mildly anisotropic, and it is perhaps reasonable to 
approximate it as completely isotropic bulk and shear viscosities with the following values: 

ri= 0.020 ±0.003 

(25) 

Cg (0.018,0.056). 

The reason for ( not being well determined by the viscosity coefficients fl2ip is that 
those coefficients do not exactly correspond to an isotropic viscosity (see equation ([2i 



5. Conclusion 



We have extended the analysis of iPenev et al.l (120071 ) to calculate effective viscosities in 
the surface convective zones of three main sequence stars: O.775M0, O.85M0 and the present 
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day Sun. We have also modified the calculation to properly account for all normalization 
factors. 

The effective viscosity we find (given by equations (120|) and (l2Ti) ) scales linearly with 
the period of the external perturbation, wit h the shear viscosity being smaller than the 



linear scaling proposed by 



ZahnI (Il966 . 



19891 ) by a factor of about ten, but in addition there 



is a significant bulk viscosity, which is assumed zero in Zahn's prescription. 

This factor in practice do es not have a dra matic effect on the tidal circularization 



period, which scales as K^^^^ ( IZahnlll966l . Il989l ). So, assuming that the above effective 



viscosity is correct in the range of per iods a p plicab 



saturation period is 2r as assumed by 



ZahnI (11966 



e to stellar binary orbits, and that the 



19891 ). the circularization cut-off period 



based on our viscosity would be within about 30% of the prediction with Zahn's scaling. 

The important difference between this effective viscosity and equation ([1]) is the 
presence of a significant bulk viscosity, the possibility that the effective viscosity is not 
isotropic, and that the linear scaling should apply only for a limited range of frequencies. 

The applicability of this result is limited by two factors: the range of applicability of 
the perturbative expansion (Eq. [7j) and the limits of the numerical simulations. 

The external shear velocities are assumed, by the perturbative expansion, to be small 
compared to the typical convective flow, and the period of the external shear should be 
neither too long nor too short. 

On one hand the limited spatial resolution of the numerical simulations means that 
only sufficiently large turbulent eddies are captured, which implies that our results do 
not apply to external forcing with very short period, for which the dissipation may be 
dominated by eddies that are too small to be reliably simulated. However, sufficiently short 
periods fall within the inertial subrange where Kolmogorov scaling holds and in that case 
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the sa me perturbational calcu lation predicts quadratic scaling of the effective viscosity with 



period (iGoodman &: Oh 



19971). 



On the other hand, the perturbative expansion we use, assumes that the perturbation 
period (T) is small compared to the turnover time (r) of the largest local eddies. In 
particular, we expect that the effective viscosity should reach a maximum value for some 
perturbation period on the order of r, and remain the same for all longer periods. This 
saturation cannot be captured by our perturbative approach since it is due to the neglected 
higher order terms. 



Penev et al. 



(120081 ) used a spectral, anelastic, ideal gas convective box, which includes 
the external forcing as part of the equations of motion, to find the effective viscosity directly 
without a perturbative treatment. They confirm that the slope of the perturbative viscosity 
is consistent with the directly calculated values in the range of its applicability, although 
for the X — z component they observe a period independent offset between the perturbative 
and direct viscosity which acts to increase the anisotropy. They also find linear scaling of 
the effective viscosity with period that saturates for T > 2tc. The magnitude of the effective 
viscosity they find based on the perturbative calculation described above is approximately 
a factor of t wo larger than the r esults presented in this paper. However, this is due to the 



fact that the 



Penev et al 



(120081) convective zone has a mixing length parameter of about 
3: double the value usually assumed for the Sun and appropriate for the simulations used 
above. 

We would like to thank the anonymous referee for detailed discussion of the 
parametrization of the viscosity which improved the paper considerably. We would also like 
to acknowledge much helpful advice that generally improved the quality of this work from 
Dr. Jeremy Goodman. 
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Table 1: The physical characteristics of the three simulations used to derive effective viscosi- 
ties. The units of Tes are K and the units of g are cms~^. 



Model Mass 


IM© 


0.85 Mq 


0.775 Mq 


Age (Gyr) 


4.55 


7 


2 


log Teff 


3.761 


3.685 


3.708 


log g 


4.44 


4.592 


4.592 


Size (Mm) {L^ x Ly x L^) 


5.4^ X 2.8 


2.72 X 1.8 


2.92 X 1.9 


Grid {N^ x Ny x N^) 


1142 X 170 


115^ X 170 


1152 X 170 


R/R© 


1.0 


0.737 


0.772 
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